home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
fft.lha
/
fft
/
vfft.cc
< prev
next >
Wrap
C/C++ Source or Header
|
1993-08-08
|
13KB
|
472 lines
// This may look like C code, but it is really -*- C++ -*-
/*
************************************************************************
*
* Verify the Fast Fourier Transform Package
*
************************************************************************
*/
#include "fft.h"
#include <math.h>
#include <builtin.h>
#include <ostream.h>
/*
*------------------------------------------------------------------------
* Timing the program execution
*/
#include <time.h>
static tms clock_acc;
static void start_timing(void)
{
times(&clock_acc);
}
static void print_timing(const char * header)
{
register int old_tick = clock_acc.tms_utime;
register double timing = (times(&clock_acc),clock_acc.tms_utime
- old_tick )/60.; // In secs
printf("\nIt took %.2f sec to perform %s\n",timing,header);
}
/*
*-----------------------------------------------------------------------
*/
// Simplified printing a vector
static void print_seq(const char * header, const Vector& v)
{
register int i;
cout << "\n" << header << "\t";
for(i=v.q_lwb(); i<=v.q_upb(); i++)
cout << form("%7.4f ",v(i));
cout << "\n";
}
static void verify_vector_identity(const Vector& v1, const Vector& v2)
{
register imax = 0;
register double max_dev = 0;
register int i;
are_compatible(v1,v2);
for(i=v1.q_lwb(); i<=v1.q_upb(); i++)
{
register double dev = abs(v1(i)-v2(i));
if( dev >= max_dev )
imax = i, max_dev = dev;
}
if( max_dev == 0 )
return;
if( max_dev < 1e-5 )
message("Two #%d elements of the vectors with values %g and %g\n"
"differ the most, though the deviation %g is small\n",
imax,v1(imax),v2(imax),max_dev);
else
_error("A significant difference between the vectors encountered\n"
"at #%d element, with values %g and %g",
imax,v1(imax),v2(imax));
}
/*
*-----------------------------------------------------------------------
* Check FFT of the Arithmetical Progression sequence
* x[j] = j
* Analytical transform is
* SUM{ j*W^(kj) } = N/(W^k - 1), k > 0,
* N*(N-1)/2, k = 0
*/
void test_ap_series(const int N)
{
cout << "\n\nVerify the computed FFT of the AP series x[j]=j\n";
cout << "j = 0.." << N-1 << "\n";
Vector xre(0,N-1);
Vector xim(xre);
register int j;
for(j=0; j<N; j++)
xre(j) = j;
xim = 0;
Vector xfe_re(xre), xfe_im(xre), xfe_abs(xre); // Exact transform
xfe_re(0) = xfe_abs(0) = N*(N-1)/2;
xfe_im(0) = 0;
register int k;
for(k=1; k<N; k++)
{
Complex arg(0,-2*PI/N * k);
Complex t = N / ( exp(arg) - 1 );
xfe_re(k) = t.real();
xfe_im(k) = t.imag();
xfe_abs(k) = abs(t);
}
FFT fft(N);
Vector xf_re(xre), xf_im(xre), xf_abs(xre);
cout << "\nPerforming Complex FFT of AP series (IM part being set to 0)\n";
start_timing();
fft.input(xre,xim);
fft.real(xf_re);
fft.imag(xf_im);
print_timing("Complex Fourier transform");
fft.abs(xf_abs);
cout << "Verifying the Re part of the transform ...\n";
verify_vector_identity(xfe_re,xf_re);
cout << "Verifying the Im part of the transform ...\n";
verify_vector_identity(xfe_im,xf_im);
cout << "Verifying the power spectrum ...\n";
verify_vector_identity(xfe_abs,xf_abs);
Vector xfr_re(xre), xfr_im(xre), xfr_abs(xre);
cout << "\nPerforming FFT of a REAL AP sequence\n";
start_timing();
fft.input(xre);
fft.real(xfr_re);
fft.imag(xfr_im);
print_timing("\"Real\" Fourier transform");
fft.abs(xfr_abs);
cout << "Check out that \"Real\" and Complex FFT give identical results";
verify_vector_identity(xfr_re,xf_re);
verify_vector_identity(xfr_im,xf_im);
verify_vector_identity(xfr_abs,xf_abs);
cout << "\nDone\n";
}
/*
*-----------------------------------------------------------------------
* Check out the orthogonality of the basis functions of FFT
*
* x[j] = W^(-l*j)
* SUM{ x[j] * W^(kj) } = 0, k <> l
* N, k = l
*/
void test_orth(const int N, const int l)
{
cout << "\n\nVerify the computed FFT for x[j] = W^(-l*j)\n";
cout << "j = 0.." << N-1 << ", l=" << l << "\n";
Vector xre(0,N-1);
Vector xim(xre);
register int j;
for(j=0; j<N; j++)
{
Complex arg(0, 2*PI/N * l * j);
Complex t = exp(arg);
xre(j) = t.real(), xim(j) = t.imag();
}
Vector xfe_re(xre), xfe_im(xre); // Exact transform
xfe_re(l) = N; xfe_im = 0;
FFT fft(N);
Vector xf_re(xre), xf_im(xre), xf_abs(xre);
cout << "\nPerforming Complex FFT\n";
start_timing();
fft.input(xre,xim);
fft.real(xf_re);
fft.imag(xf_im);
print_timing("Complex Fourier transform");
fft.abs(xf_abs);
cout << "Verifying the Re part of the transform ...\n";
verify_vector_identity(xfe_re,xf_re);
cout << "Verifying the Im part of the transform ...\n";
verify_vector_identity(xfe_im,xf_im);
cout << "Verifying the power spectrum ...\n";
verify_vector_identity(xfe_re,xf_abs);
cout << "\nDone\n";
}
/*
*-----------------------------------------------------------------------
* Check FFT of the truncated Arithmetical Progression sequence
* x[j] = j, j=0..N/2
* Analytical transform is
* SUM{ j*W^(kj) } = N/2 (W^k - 1), k > 0 and even
* 2*W^k/(W^k - 1)^2 - N/2 * 1/(W^k - 1), k being odd
* N/2 * (N/2-1)/2, k = 0
*/
void test_ap_series_pad0(const int N)
{
cout << "\n\nVerify the computed FFT of the truncated AP sequence x[j]=j\n";
cout << "j = 0.." << N/2-1 << ", with N=" << N << "\n";
Vector xre(0,N/2-1);
Vector xim(xre);
register int j;
for(j=0; j<N/2; j++)
xre(j) = j;
xim = 0;
Vector xfe_re(0,N-1), xfe_im(0,N-1), xfe_abs(0,N-1); // Exact transform
xfe_re(0) = xfe_abs(0) = N/2 * (N/2 - 1)/2;
xfe_im(0) = 0;
register int k;
for(k=1; k<N; k++)
{
Complex wk = exp(Complex(0,-2*PI/N * k));
Complex t = 1 / ( wk - 1 );
if( k & 1 )
t = 2*wk*t*t - N/2 * t; // for k odd
else
t *= N/2; // for k even
xfe_re(k) = t.real();
xfe_im(k) = t.imag();
xfe_abs(k) = abs(t);
}
FFT fft(N);
Vector xf_re(0,N-1), xf_im(0,N-1), xf_abs(0,N-1);
cout << "\nPerforming Complex FFT (with IM part being set to 0)\n";
start_timing();
fft.input_pad0(xre,xim);
fft.real(xf_re);
fft.imag(xf_im);
print_timing("Complex Fourier transform");
fft.abs(xf_abs);
if( N <= 16 )
{
print_seq("Source Vector ",xre);
print_seq("Computed cos transform",xf_re);
print_seq("Computed sin transform",xf_im);
}
cout << "Verifying the Re part of the transform ...\n";
verify_vector_identity(xfe_re,xf_re);
cout << "Verifying the Im part of the transform ...\n";
verify_vector_identity(xfe_im,xf_im);
cout << "Verifying the power spectrum ...\n";
verify_vector_identity(xfe_abs,xf_abs);
Vector xfr_re(xf_re), xfr_im(xf_re), xfr_abs(xf_re);
cout << "\nPerforming FFT of a REAL AP sequence\n";
start_timing();
fft.input_pad0(xre);
fft.real(xfr_re);
fft.imag(xfr_im);
print_timing("\"Real\" Fourier transform");
fft.abs(xfr_abs);
cout << "Check out that \"Real\" and Complex FFT give identical results\n";
verify_vector_identity(xfr_re,xf_re);
verify_vector_identity(xfr_im,xf_im);
verify_vector_identity(xfr_abs,xf_abs);
Vector xfh_re(xre), xfh_im(xre), xfh_abs(xre);
cout << "Check out the functions returning the half of the transform\n";
fft.real_half(xfh_re);
fft.imag_half(xfh_im);
fft.abs_half(xfh_abs);
Vector xfhr_re(xre), xfhr_im(xre), xfhr_abs(xre);
for(j=0; j<N/2; j++)
xfhr_re(j) = xfh_re(j), xfhr_im(j) = xfh_im(j), xfhr_abs(j) = xfh_abs(j);
verify_vector_identity(xfhr_re,xfh_re);
verify_vector_identity(xfhr_im,xfh_im);
verify_vector_identity(xfhr_abs,xfh_abs);
cout << "\nDone\n";
}
/*
*-----------------------------------------------------------------------
* Verify the sin/cos Fourier transform
*
* r*exp( -r/a ) <=== sin-transform ===> 2a^3 k / (1 + (ak)^2)^2
* r*exp( -r/a ) <=== cos-transform ===> a^2 (1 - (ak)^2)/(1 + (ak)^2)^2
*/
static void test_lorentzian(void)
{
const double R = 20; // Cutoff distance
const int n = 512; // No. of grids
const double a = 4; // Constant a, see above
const double dr = R/n, // Grid meshes
dk = PI/R;
cout << "\n\nVerify the sin/cos transform for the following example\n";
cout << "\tr*exp( -r/a )\t<=== sin-transform ===>\t2a^3 k/(1 + (ak)^2)^2\n";
cout << "\tr*exp( -r/a )\t<=== cos-transform ===>\t"
"a^2 (1-(ak)^2)/(1 + (ak)^2)^2\n";
cout << "\nParameter a is " << form("%.2f",a);
cout << "\nNo. of grids " << n;
cout << "\nGrid mesh in the r-space dr = " << form("%.3f",dr);
cout << "\nGrid mesh in the k-space dk = " << form("%.3f",dk);
FFT fft(2*n,dr);
cout << "\nCheck out the inquires to FFT package about N, dr, dk, cutoffs\n";
assert( fft.q_N() == 2*n );
assert( fft.q_dr() == dr );
assert( fft.q_dk() == dk );
assert( fft.q_r_cutoff() == R );
assert( fft.q_k_cutoff() == dk*n );
Vector xr(0,n-1); // Tabulate the source function
register int j;
for(j=0; j<n; j++)
{
double r = j*dr;
xr(j) = r * exp( -r/a );
}
Vector xs(xr), xc(xr);
start_timing();
fft.sin_transform(xs,xr);
print_timing("Sine transform");
start_timing();
fft.cos_transform(xc,xr);
print_timing("Cosine transform");
Vector xs_ex(xs), xc_ex(xc); // Compute exact transforms
for(j=0; j<n; j++)
{
double k = j * dk;
Complex t = 1/Complex(1/a,-k);
t = t*t;
xc_ex(j) = t.real();
xs_ex(j) = t.imag();
}
compare(xs,xs_ex,"Computed and Exact sin-transform");
compare(xc,xc_ex,"Computed and Exact cos-transform");
Vector xt(xc); xt = xc; xt -= xc(xc.q_upb());
compare(xt,xc_ex,"Computed cos-transform with DC component removed, and "
"exact result");
Vector xs_inv(xr), xc_inv(xr); // Compute Inverse transforms
start_timing();
fft.sin_inv_transform(xs_inv,xs);
print_timing("Inverse sine transform");
start_timing();
fft.cos_inv_transform(xc_inv,xc);
print_timing("Inverse cosine transform");
compare(xs_inv,xr,"Computed inverse sin-transform vs the original function");
compare(xc_inv,xr,"Computed inverse cos-transform vs the original function");
xt = xc_inv; xt -= xc_inv(xc_inv.q_upb());
compare(xt,xr,"Computed inverse cos-transform with DC component removed,\n"
"and the original function");
cout << "\nDone\n";
}
/*
*-----------------------------------------------------------------------
* Verify the cos Fourier transform
*
* exp( -r^2/4a ) <=== cos-transform ===> sqrt(a*pi) exp(-a*k^2)
*/
static void test_gaussian(void)
{
const double R = 20; // Cutoff distance
const int n = 512; // No. of grids
const double a = 4; // Constant a, see above
const double dr = R/n, // Grid meshes
dk = PI/R;
cout << "\n\nVerify the sin/cos transform for the following example\n";
cout << "\texp( -r^2/4a )\t<=== cos-transform ===> sqrt(a*pi) exp(-a*k^2)\n";
cout << "\nParameter a is " << form("%.2f",a);
cout << "\nNo. of grids " << n;
cout << "\nGrid mesh in the r-space dr = " << form("%.3f",dr);
cout << "\nGrid mesh in the k-space dk = " << form("%.3f",dk);
FFT fft(2*n,dr);
cout << "\nCheck out the inquires to FFT package about N, dr, dk, cutoffs\n";
assert( fft.q_N() == 2*n );
assert( fft.q_dr() == dr );
assert( fft.q_dk() == dk );
assert( fft.q_r_cutoff() == R );
assert( fft.q_k_cutoff() == dk*n );
Vector xr(0,n-1); // Tabulate the source function
register int j;
for(j=0; j<n; j++)
{
double r = j*dr;
xr(j) = exp(-r*r/4/a );
}
Vector xc(xr);
start_timing();
fft.cos_transform(xc,xr);
print_timing("Cosine transform");
Vector xc_ex(xr); // Compute exact transforms
for(j=0; j<n; j++)
{
double k = j * dk;
xc_ex(j) = sqrt(PI*a) * exp( -a*k*k );
}
compare(xc,xc_ex,"Computed and Exact cos-transform");
xc -= xc(xc.q_upb());
compare(xc,xc_ex,"Computed with DC removed, and Exact cos-transform");
Vector xc_inv(xr); // Compute Inverse transforms
start_timing();
fft.cos_inv_transform(xc_inv,xc);
print_timing("Inverse cosine transform");
compare(xc_inv,xr,"Computed inverse cos-transform vs the original function");
xc_inv -= xc_inv(xc_inv.q_upb());
compare(xc_inv,xr,"Computed inverse with DC removed vs the original");
cout << "\nDone\n";
}
/*
*------------------------------------------------------------------------
* Root module
*/
main()
{
cout << "\n\n" << _Equals <<
"\n\n\t\tVerify Fast Fourier Transform Package\n\n";
test_ap_series(8);
test_ap_series(1024);
test_orth(1024,1);
test_ap_series_pad0(16);
test_ap_series_pad0(1024);
test_lorentzian();
test_gaussian();
}